function x = custom_finv(p, v1, v2, pin)
% custom_finv  - returns inverse of the F-distribution
%
% FORMAT:       x = custom_finv(p, v1, v2);
%
% Input fields:
%
%       p           CDF of F distribution with v1, v2 d.f. at points x
%       v1, v2      d.f. of distribution function
%
% Output fields:
%
%       x           F-variate

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 3 || ...
   ~isnumeric(p) || ...
   ~isnumeric(v1) || ...
   ~isnumeric(v2) || ...
    isempty(v1) || ...
    isempty(v2) || ...
    any(isnan(p(:))) || ...
    any(isinf(v1(:)) | isnan(v1(:)) | v1(:) < 0) || ...
    any(isinf(v2(:)) | isnan(v2(:)) | v2(:) < 0)
    error( ...
        'BVQXtools:BadArgument', ...
        'Missing or invalid argument given.' ...
    );
end
if isempty(p)
    x = [];
    return;
end
osize = size(p);
if numel(p) == 1
    if numel(v1) > 1
        osize = size(v1);
    elseif numel(v2) > 1
        osize = size(v2);
    end
end
if ~isa(p, 'double')
    p = double(p(:));
else
    p = p(:);
end
if ~isa(v1, 'double')
    v1 = double(v1(:));
else
    v1 = v1(:);
end
if ~isa(v2, 'double')
    v2 = double(v2(:));
else
    v2 = v2(:);
end

% create output
x = zeros(osize);
x(p < 0 | p > 1 | isnan(p) | ~(v1 > 0) | ~(v2 > 0)) = NaN;
x(p == 1 & v1 > 0 & v2 > 0) = Inf;

% regular cases
k = find(p > 0 & p < 1 & v1 > 0 & v1 > 0);
if ~isempty(k)
    if nargin > 3 && ...
        islogical(pin) && ...
        numel(pin) == 1 && ...
        pin
        if numel(v1) == 1 && ...
            numel(v2) == 1
            x(k) = ((1 ./ custom_betainv(p(k), v2 / 2, v1 / 2) - 1) * v2 ./ v1);
        elseif numel(v1) == 1
            x(k) = ((1 ./ custom_betainv(p(k), v2(k) ./ 2, v1 / 2) - 1) .* ...
                v2(k) ./ v1);
        elseif numel(v2) == 1
            x(k) = ((1 ./ custom_betainv(p(k), v2 / 2, v1(k) ./ 2) - 1) .* ...
                v2 ./ v1(k));
        else
            x(k) = ((1 ./ custom_betainv(p(k), v2(k) ./ 2, v1(k) ./ 2) - 1) .* ...
                v2(k) ./ v1(k));
        end
    else
        if numel(v1) == 1 && ...
            numel(v2) == 1
            x(k) = ((1 ./ custom_betainv(1 - p(k), v2 / 2, v1 / 2) - 1) * v2 ./ v1);
        elseif numel(v1) == 1
            x(k) = ((1 ./ custom_betainv(1 - p(k), v2(k) ./ 2, v1 / 2) - 1) .* ...
                v2(k) ./ v1);
        elseif numel(v2) == 1
            x(k) = ((1 ./ custom_betainv(1 - p(k), v2 / 2, v1(k) ./ 2) - 1) .* ...
                v2 ./ v1(k));
        else
            x(k) = ((1 ./ custom_betainv(1 - p(k), v2(k) ./ 2, v1(k) ./ 2) - 1) .* ...
                v2(k) ./ v1(k));
        end
    end
end

% reshape
x = reshape(x, osize);
